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ABSTRACT 

The efficacy of the Hill stability (HS) criterion is compared to other known chaos 
indicators such as the maximum Lyapunov exponent (MLE) and Mean Exponential 
Growth factor of Nearby Orbits (MEGNO) maps. The orbits of four individual planets 
in four known binary star systems, 7 Cephei, Gliese 86, HD 41004, and HD 196885, 
are numerically integrated using various numerical techniques to assess the chaotic 
or quasi-periodic nature of the dynamical system considered. The Hill stability which 
measures the orbital perturbation of a planet around the primary star due to the 
secondary star is calculated for each system. The maximum Lyapunov exponent time 
series are generated to measure the divergence/convergence rate of stable manifolds, 
which are used to differentiate between chaotic and non-chaotic orbits. Then, we cal- 
culate dynamical MEGNO maps from solving the variational equations along with the 
equations of motion. These maps allow us to accurately differentiate between stable 
and unstable dynamical systems. The results obtained from the analysis of HS, MLE, 
and MEGNO maps are analysed for their dynamical variations and resemblance. The 
qualitative efficiency of each indicator is analysed which demonstrates that HS can be 
used as an alternative to MLE. The HS test for the planets shows stability and quasi- 
periodicity for at least ten million years. The MLE and the MEGNO maps have also 
indicated the local quasi- periodicity and global stability in relatively short integration 
period. Based on our discussion, the HS criterion is found to be a comparably efficient 
tool to measure the stability of planetary orbits with respect to different simulation 
timespans. 
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1 INTRODUCTION 

The discovery of extra solar planets has been growing sub- 
stantially since the first planet, 51 Pegasi b, was detected 
almost two decades ago (Mayor & Queloz 1995). Since then 
843 extra solar planets have been confirmed as of Octo- 
ber 22, 2012^. Near half of solar type stars (Duquennoy 
& Mayor 1991; Raghavan et al. 2006) and a third of all 
stars in the Galaxy (Raghavan et al. 2010) are in binary 
or multi-star system with 40 planets confirmed in such sys- 
tems (Desidera & Barbieri 2007). The confirmation of the 
existence of planets in binaries has raised a new astrophys- 
ical challenge which includes the study of long term orbital 
stability of such planets. The ultimate destiny of exoplanet 
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research, including observations from the Kepler mission^, 
is to detect a planet on a stable orbit within the habitable 
zone of the host star (Borucki et al. 1997; Koch et al. 2007; 
Borucki et al. 2008, 2010). The degree of stability is largely 
governed by the planet's semi-major axis, eccentricity and 
orbital inclination. Orbital long-term stability is believed to 
be a necessary condition for life to develop. Our aim is to 
find an answer to the century old question, "Are we alone 
in the Milky Way Galaxy?". 

By using different stability criteria the question of sta- 
bility has been addressed by many others in the past. While 
studying the Trojan type orbits around Neptune, Zhou et al. 
(2009) showed that the inclination of orbits can be as high 
as 60° while maintaining orbital stability. Several authors 
(Szebehely 1980; Szebehely & McKenzie 1981; Szenkovits 
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& Mako 2008) calculated orbital stability of planets by us- 
ing several techniques that include the integrals of motion, 
zero velocity surfaces (ZVS), and a Hill stable region that 
is mapped by a parameter space of orbital radius and mass 
ratio, fi, for a coplanar circular restricted three body prob- 
lem (CRTBP). Quarles et al. (2011) has also used the max- 
imum Lyapunov exponent (MLE) to determine the orbital 
stability or instability for the CRTBP case. The stability 
limits were defined based on the values of MLE that are de- 
pendent on the mass ratio /i of the binaries and the initial 
distance ratio po of the planet. Other chaos indicator tech- 
niques such as Mean Exponential Growth factor of Nearby 
Orbits (MEGNO) maps have also been used to study the 
dynamical stability of irregular satellites (Hinse et al. 2010) 
and extrasolar planet dynamics (Gozdziewski & Maciejew- 
ski 2001; Gozdziewski et al. 2001). The MEGNO criterion is 
known to be efficient in distinguishing between chaotic and 
quasi-periodic initial conditions within a dynamical phase 
space. 

Knowing the orbital stability of planets is a crucial step 
for further studies of planetary systems. In order to probe 
a planet for its habitability, there exists a primary require- 
ment that the system be orbitally stable. In this paper, we 
have used three well known chaos indicators HS, MLE, and 
MEGNO maps in the study of the orbital perturbation of 
planets in the selected stellar binaries. We have compared 
all the results from three different methods. The compari- 
son study allows us to determine the best tool for a stability 
analysis of the considered systems. 

This paper is outlined as follows. In Section 2 we discuss 
the basic theory of the analysis tools. In Section 3 we present 
our numerical methods and the results from the comparison 
of chaos indicators followed by discussion. Finally, we will 
conclude in Section 4 with a brief overview of our results. 



2 THEORY 

2.1 Basic Definitions and Equations 

For the motion of a planet of mass mi around a star of mass 
777-2 in an orbital plane, the initial position of the planet at 
any given time is given by (Murray & Dermott 1999) 

(X\ / cos Q cos (cj + /) — sin Q cos i sin {oj + /)\ 

Y \ = r sin O cos (cj + /) + cos Q cos i sin (cj + /) , 
Z J \ sin i sin (a; + /) / 

(1) 

where the orbital parameters (r, z, cj, Q, /) are the radius 
vector, inclination angle, argument of periapsis, longitude of 
ascending node and true anomaly, respectively. 

Then we can calculate the velocities in the x, y and z 
directions by taking the time derivative of Eq.(l) and obtain 

Vx = A{e sin / [cos Q cos (cj + /) — sin Q cos i sin (cj + /)] 
— (1 + e cos /) [cos Q sin (cj + /) + sin Q cos i cos (cj + /)]}, 

Vy = A{e sin / [sin Q cos (cj + /) + cos Q cos i sin {oj + /)] 
+ (1 + e cos /) [cos Q cos i cos (cj + /) — sin Q sin (cj + /)]}, 

Vz = Asini [cos (u + /)(1 + ecos/) + sin (O + /)esin/] . 

(2) 



Equations (2) have been simplified using the following 
relations 



f- 

A: 



rVl — 

na 



(1 + ecos/), 



(3) 



where (/, n) are the time derivative of true anomaly 
and mean motion respectively. For our special case we con- 
sidered the inclination i = and while calculating the initial 
conditions we used the values of other parameters (Q, oj and 
/) whenever they have been observationally determined. 

Our particular interest is on the stellar binaries that 
are less than or equal to 25 AU apart. For the stars with 
greater than 25 AU separations, the effects of the secondary 
star on the planet would not be significant, especially while 
considering the intent of our present study. 

The list of planets in the binaries and their orbital pa- 
rameters are given in Tables 1 through 4. The initial setup in 
our simulations is in a barycentric coordinate system with 
Star B to the left and Star A to the right. The positive 
X-axis is taken to be the reference. The true anomaly and 
Right Ascension of ascending node of both the binary and 
planet are assumed to be equal to zero. 

For two primaries in elliptic orbits moving about their 
barycentre the dynamical system of a third smaller mass 
can be written as the first order differential equations. 
The equations of motion are given below (Szebehely 1967; 
Szenkovits & Mako 2008) 



X = u u = 2v -\- 



(1 -h ecos /) 



a(x -h /i) fi(x — 1 -h /i) 



y = V V = —2u -h 
z = w w = —z -\- 
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(1 -h ecos /) 
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(777-1 + ^77-2) ' 

a = 1 — fi, 

2 / n2 . 2 , 2 

-1 = -\-y + z , 

2 / I n2 , 2 I 2 

-2 = [x + a) +y + z . 



(5) 



The Jacobi constant (Co) for a initial conditions {xo, 
Vo, Zo, fo) is given by 



Co 



1 -h ecos(/o) 



. 2 



.2 

Vo 



.2 
Zq. 



(6) 



In equations 4, the variables represent the velocity of 
a test particle (planet) in Cartesian coordinates {x, y, z). 
The distances ri and r2 are defined in terms of mass ratio, 
normalised coordinates, and the position of the stars within 
a rotating coordinate system. 

2.2 Lyapunov Exponents 

For stable planetary orbits, the two nearby trajectories 
in phase space will converge and for unstable orbits, the 
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trajectories diverge exponentially. The rate of divergence 
is measured by using the method of Lyapunov exponents 
(Lyapunov 1907). Wolf et al. (1985) developed a numeri- 
cal method of computing the Lyapunov exponents in FOR- 
TRAN following the earher works by Benettin et al. (1980). 

Lyapunov exponents are commonly used because they 
give the measure of an attractor of a dynamical system as it 
converges or diverges in phase space. The positive Lyapunov 
exponents measure the rate of divergence of neighboring or- 
bits, whereas negative exponents measure the convergence 
rates between stable manifolds (Tsonis 1992; Ott 1993). The 
sum of all Lyapunov exponents is less than zero for dissipa- 
tive systems (Musielak & Musielak 2009) and zero for non- 
dissipative (Hamiltonian) systems (Hilborn & Sprott 1994). 
Lyapunov exponents for the circular restricted three body 
problems (CRTBP) have been calculated previously (Gonczi 
& Froeschle 1981; Murray & Holman 2001). In this work, we 
have calculated the Lyapunov exponents for the elliptic re- 
stricted three body problem (ERTBP) and one case of the 
elliptic restricted four body problem (HD 41004). In order to 
calculate the Lyapunov exponents a dynamical system with 
n degrees of freedom is represented in a 2n phase space. 
Then for the case of ERTBP, state vectors {2n) containing 
6 elements are used to calculate the Lyapunov exponents. 
The details on the calculation of Jacobian J from the equa- 
tions of motion can be found in Quarles et al. (2011). 

For a Hamiltonian system (see above) to be stable, the 
sum of all the Lyapunov exponents should be zero. In or- 
der to numerically meet such a criterion, a simulation of the 
system would require an impractically long period of time. 
Within the limits of simulation, the sum of all six exponents 
must remain numerically close (~ 10~^°)to zero. We have 
used the largest positive Lyapunov exponent to determine 
the magnitude of the chaos while ignoring the lesser positive 
and the negative exponents as they do not provide signifi- 
cant additional information about the evolving system. The 
positive maximum Lyapunov exponent is known to indicate 
a chaotic behaviour in both dissipative (Hilborn & Sprott 
1994) and non-dissipative (Ozorio de Almeida 1990) sys- 
tems. Chaos can be proven up to the integration time if a 
given chaos indicator has converged to an unstable mani- 
fold. However, quasi-periodicity or regular motion can only 
be proven up to the considered integration time. Quasi- 
periodicity or stable motion can never rigorously be prooven 
for the 3 (or more) body system. We never know how the 
system might evolve after the considered integration time. 



2.3 Hill Stability 

Hill (1878a,b) developed the equations of motion for a par- 
ticle around the primary mass in presence of a nearby sec- 
ondary mass. The purpose of the Hill equations was to cal- 
culate the orbital perturbation of the particle due to the 
secondary mass. Later the idea was further developed and 
used in the study of orbital stability of planets (Szebehely 
1967; Walker & Roy 1981; Marchal & Bozis 1982). 

The significant radial gravitational influence of the sec- 
ondary mass reaches as far as the Lagrange points, LI and 
L2, forming the Hill sphere (Hill 1878a, b). The contour lines 
within the sphere are the zero velocity curves. After mea- 
suring a particle's position and velocity a constant of mo- 
tion relation can be implemented (Szebehely 1967; Murray 



7 Cephei 


A 


B 


Ab 


Mass 

Semimajor Axis (a) 

Eccentricity (e) 

Argument of Periapsis {uop) 


1.4 M0 
19.02 AU 
0.4085 
0° 


0.362 M( 

180° 


© 1.6 Mj 

1.94 AU 
0.115 

94° 


Table 1. Orbital parameters of 7 Cephei (Neuhauser et. al. 2007) 


Gliese 86 


A 


B 


Ab 


Mass 

Semimajor Axis (a) 

Eccentricity (e) 

Argument of Periapsis {uOp) 


0.79 M0 
18.42 AU 
0.3974 
0° 


0.55 M© 
180° 


3.91 Mj 

0.113 AU 
0.0416 
269° 


Table 2. Orbital parameters of Gliese 86 (Lagrange 
Butler et al. 2006; Santos et al. 2004) 


et al. 2006; 


HD 196885 


A 


B 


Ab 


Mass 

Semimajor Axis (a) 

Eccentricity (e) 

Argument of Periapsis {ujp) 


1.33 M© 
21 AU 
0.42 
0° 


0.45 M© 
180° 


2.98 Mj 
2.6 AU 
0.48 
93.2° 


Table 3. Orbital parameters of HD 196885, (Chauvin 


et al. 2007) 


HD 41004 


B 


Ab 


Bb 


Mass 

Semimajor Axis (a) 

Eccentricity (e) 

Argument of Periapsis {00 p) 


0.4 M© 
22 AU 
0.065 
180° 


2.3 Mj 
1.31 AU 
0.39 

97° 


18.64 Mj 
0.0177 AU 
0.081 
178.5° 



Table 4. Orbital parameters of HD 41004, M^ = 0.7M© (Zucker 
et al. 2004; Butler et al. 2006) 



& Dermott 1999) 2U - v'^ ^ Cj, where v is the velocity, U 
is the generalised potential, and Cj is the constant of in- 
tegration called the Jacobi constant. When the velocity of 
the particle is zero, 2U = Cj, a contour represents a zero 
velocity surface (ZVS) and the motion of a particle within 
such a surface is considered Hill stable. 

The measure of Hill stability, *S'(/), is given by a pa- 
rameter dependent potential, Q(x, /), where / is the true 
anomaly and Ccr is the value of the Jacobi constant at the 
Hill radius or the Lagrange point LI (Szenkovits & Mako 
2008) 



■ ez' cos(/)] + — 



+ 



M), 

[2Q(x,/). 



[l + ecos(/)]^'] -1. (7) 



Using the orbital parameters obtained from the numeri- 
cal integration the potential, Q(x, /), is calculated to obtain 
the Hill stability function 5'(/). Although the Hill stability 
function depends on the true anomaly, it can also be repre- 
sented as a time series. We have implemented this represen- 
tation in our results concerning the Hill stability function. 
When the measure of *S'(/) of a planet is positive then we 
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have the indication of quasi-periodic orbits and its motion 
is Hill stable. But when the measure of S{f) is negative 
then the planet enters the chaotic region, hence loosing its 
stability. 

2.4 The MEGNO Chaos Indicator 

The MEGNO criterion was first introduced by Cincotta & 
Simo (1999, 2000); Cincotta et al. (2003) and found wide- 
spread applications in dynamical astronomy (Gozdziewski 
et al. 2001, 2008; Gozdziewski & Migaszewski 2009; Hinse 
et al. 2008, 2010; Frouard et al. 2011; Compere et al. 2012; 
Rostov et al. 2012). The MEGNO (usually denoted as {¥)) 
formalism has the following mathematical properties. In gen- 
eral, MEGNO has the parameterisation (Y) = axt-\- 13 (see 
references above). For a quasi-periodic initial condition, we 
have a ~ 0.0 and /3 ~ 2.0 (or (Y) 2.0) for t ^ oo asymp- 
totically. If the orbit is chaotic, then (Y) — At/2 for t — cxd. 
Here A is the maximum Lyapunov exponent (MLE) of the or- 
bit. In practice, when generating our MEGNO maps, we ter- 
minate a given numerical integration of a chaotic orbit when 
(Y) > 12.0. Quasi-periodic orbits have |(y) - 2.0| = 0.001. 

We used the MECHANIC^ software (Slonina et al. 
2012) to calculate the MEGNO maps on a multi-CPU com- 
puting environment. Typically we allocated 60 CPUs for 
the calculation of one map considering a typical grid of 
(500 X 300) initial conditions in (a, e) space. The numeri- 
cal integration of the equations of motion and the associ- 
ated variational equations (Mikkola & Innanen 1999) are 
based on the ODEX integration software (Hairer et al. 
1993) which implements a Gragg-Bulirsch-Stoer algorithm. 
The MEGNO indicator is calculated from solving two addi- 
tional differential equations as outlined in Gozdziewski et al. 
(2001). 

MEGNO and the maximum Lyapunov exponent (MLE) 
have a close relation and provide the magnitude of the ex- 
ponential divergence of orbits. Froeschle et al. (1997) in- 
troduced the fast Lyapunov indicator (FLI), which exhibits 
the least dependency on initial conditions. Recently, Mestre 
et al. (2011) showed that MEGNO and FLI are related to 
each other. FLI is used to detect weak chaos and is consid- 
ered a faster means to determine the same characteristics as 
MLE. Very recently, Mafhone et al. (2011) compared various 
chaos indicators including FLI and MEGNO. The MEGNO 
technique and FLI are considered to be in the same class of 
chaos detection tools (Morbidelli 2002), and we have chosen 
the MEGNO technique to compare it versus the Hill stabil- 
ity criterion. More on mathematical properties of MEGNO 
and its relationship with the Lyapunov exponents can be 
found in Hinse et al. (2010). 



3 RESULTS AND DISCUSSION 
3.1 Numerical Simulation 

To establish the Hill stability (HS) criterion and calculate 
the maximum Lyapunov exponents (MLE), we numerically 
simulated each of the planets in the stellar binaries using a 
Yoshida sixth order symplectic and a Gragg-Bulirsch-Stoer 

^ http://www.git.astri.umk.pl/projects/mechanic 



integration scheme (Yoshida 1990; Grazier et al. 1996; Hairer 
et al. 1993). A stepping of e=10~^ years/step was used in 
each case to have a better measure of the precision of the 
integration scheme. The error in energy was calculated at 
each step which falls in the range of 10~^^ to 10~^° dur- 
ing the total integration period. Numerical simulations were 
completed for a million years to calculate the MLE and 10 
million years to calculate the HS. MEGNO maps are calcu- 
lated using 150,000 years per initial condition. 

Calculations involving chaos indicators are unaffected 
by longer simulation time. Thus, the purpose of simulating 
for 10 million years was to display the evolution of eccen- 
tricity and semi-major axis without applying any indicator 
tools to provide a consistent comparison in order to analyse 
the orbital behaviour and to establish the full effectiveness 
of the indicators. The time series plots for two of the selected 
systems, 7 Cephei and HD 196885, are relatively constant 
with minor oscillations for 10 million years (Fig. 1). In these 
cases we found that the eccentricity of the giant planets is 
oscillating with a constant amplitude. For example. Figs, 
la and Ic demonstrate oscillations from to 0.1 and 0.4 to 
0.5 in values of eccentricity for 7 Cephei and HD 196885, 
respectively. The amplitude of the oscillation changed with 
a different choice of initial conditions. As a result, specific 
choices can minimise the oscillation amplitude and can ren- 
der the simulation for 7 Cephei to be in closer agreement 
with previous studies by Haghighipour (2006). One such ini- 
tial condition involved the choice of eccentricity. If e = for 
the planetary orbit initially, we observed that the amplitude 
of oscillation is minimum and consistent with Haghighipour 
(2006) while the amplitude increases with larger initial e 
values. 

3.2 MLE: Indicator Analysis 

The maximum Lyapunov exponent (MLE) time series for 
the simulated planets in the stellar binaries are given in Figs. 
2 and 3. The MLE is plotted using a logarithmic scale along 
the y-axis and a linear scale along the time-axis. We obtained 
six Lyapunov exponents from our simulation among which 
three are negative and three are positive. We inspected the 
first three positive LEs and found the rate of change in mag- 
nitude of the largest value which is used for our purpose of 
establishing the stability of a system. 

In Figures 2 and 3, we have taken the maximum Lya- 
punov exponent as the primary indicator of the orbital sta- 
bility. For a given initial condition, the MLE must quickly 
drop below a cutoff value (Quarles et al. 2011) and decrease 
at a constant rate where we can determine it as stable or 
unstable. The MLE of the planets within each system show 
that the observed systems are in stable configurations, as ex- 
pected from the observations. In each case the MLE starts 
on the order of 10^ and slowly converges down by orders of 
10. The MLE for 7 Cephei, Fig. 2a, starts initially on the 
order of 1 and quickly drops down to -10 on a logarithmic 
scale. Then it slowly decreases to -13 in a million years. In 
the case of Gliese 86, Fig. 2b, the MLE decreases very slowly 
between -10 and -12. Considering this decreasing trend and 
the nature of Lyapunov exponents. Section 2.2, the results 
refiect the outcome of orbital stability for the planet. The 
MLE for the planet in Gliese 86 also oscillates every 0.3 mil- 
lion years which is indicative of a near resonance behaviour. 



© 2002 RAS, MNRAS 000, 1-?? 



Monthly Notices: MT^X2s guide for authors 5 



(a) 


— Giant Planet 
— y Cephei 











(b) 








— Giant Planet 
— y Cephei 









(c) 



-Giant Planet 
HD 196885 



& 








— Giant Planet 
— HD 196885 









Time(x 10° yrs) 



Time(x 10° yrs) 



Time(x 10° yrs) 



2 4 6 8 10 
Time(x 10^ yrs) 



Figure 1. Variation of orbital elements for the giant planet in 7 Cephei and HD 196885 simulated for 1 x 10^ yrs. 



Based on the MLE time series, Fig. 2c, the orbit of the planet 
in HD 196885 system is in a stable configuration. Similar to 
the case of 7 Cephei, the MLE values for HD 196885 are 
scattered below its major trend line of -13 displaying the 
possible perturbation caused by the secondary star. The os- 
cillations observed every 0.2 million years are due to the near 
resonance behaviour, a similar characteristic shown by the 
planet in Gliese 86. 

Because HD 41004 is a four body system, we examined 
it in two ways. First, we considered the system setup as an 
elliptic restricted four body problem with the primary star 
(A), secondary star (B) and two planets (Ab and Bb) or- 
biting each of the stars. Then, we separately generated the 
MLE for each of the planets (Figs. 3a and 3b). We found that 
for both of the planets, the MLE rapidly decreases until - 
10 (on a logarithmic scale) and then levels off. The orbital 
stability is ensured for both of the planets. Secondly, we re- 
duced the system to an elliptic restricted three body problem 
(ERTBP) to simplify the problem. Since the planet Bb is a 
brown dwarf, we moved the reduced mass of the planet and 
its host star to their respective bary centre. Then, with star 
A and star (B+Bb) in the system we calculated the orbital 
perturbation on the planet Ab due to the star-planet sys- 
tem (B+Bb). The MLE of the planet Ab, Fig. 3c, is slowly 
decreasing from -10 to -14, as we would expect. There is not 
a significant change in the MLE time series of the planet Ab 
when comparing between the 4 body and 3 body represen- 
tations; however, the MLE for the reduced system oscillates. 
In any case the MLE time series indicates stable orbits for 
the planets in HD 41004 binary system. The reduction al- 
lowed us to compare MLE time series for two different mod- 
els and to generate the MEGNO maps more easily, and was 
the chosen setup for further analysis in Sect. 3.4. 



3.3 Establishing the Hill Stability criterion 

The Hill stability time series for the planets in 7 Cephei, 
Gliese 86, HD 196885 and HD 41004 are shown in the Figs. 
4 and 5. For all the cases we have considered the stars and 
the planets to be coplanar (i = O'^). 

We found that the measure of Hill stability for 7 Cephei 
stays positive throughout the integration period (10 million 
years). At intervals of 2-3 and 7-8 million years a spike in 
the Hill stability is noted, which indicates that the planet 
shows a quasi-periodic motion every 5 million years. The 
oscillations on average are small and positive which reflect 
the condition for HS criteria established by Szenkovits & 
Mako (2008). Their calculation of Hill stability is positive 



and constantly increasing in time for a million years, which 
is consistent with our result. However, the plots are limited 
to a million years in the calculations of Szenkovits & Mako 
(2008) which makes the periodicity for long term oscillations 
for the planets to remain unclear. 

Similarly, the Hill stability time series stays positive and 
constant for the planet in the HD 196885 and Gliese 86 
binary systems (Figs. 4b and 4c). The planet in HD 196885 
displays a quasi periodic motion every 2 million years. For 
the planet in Gliese 86, the value of HS initially increases 
from to 2 million years and then oscillates every 0.3 million 
years with constant amplitude. The increment of HS value 
is consistent with Szenkovits & Mako (2008) for the first 
million years. 

At first the values of Hill stability for the planets in the 
HD 41004 binary system was calculated by considering the 
system as an elliptic restricted four body problem. The Hill 
stability time series of the planets are shown in Figs. 5a and 
5b. The HS values for the planets in HD 41004 Ab and Bb 
are oscillating with an amplitude ranging from 0-200 and 
0-100 respectively. Since the stability function S(f) remains 
positive for the entire simulation time, we can conclude that 
the system is Hill stable. The Hill stability function calcu- 
lated by Szenkovits & Mako (2008) for the same planets 
is positive but constantly increasing for HD 41004 Ab and 
constantly decreasing for HD 41004 Bb. In our results the 
stability values are increasing for the first half million years 
and remain constant thereafter. Szenkovits & Mako calcu- 
lated the HS for only a million years and our HD 41004 Ab 
result is consistent with their previous study. However, for 
HD 41004 Bb, the positive HS values are constantly decreas- 
ing which indicates that the planet my loose its stability 
beyond one million years. 

The HS time series for HD 41004 system indicates or- 
bital stability of the planets Ab and Bb due to the fact that 
the HS values remain positive. However, the HS values for 
both the planets are noisy and randomly fluctuating because 
it was represented as a four body system. This might be 
due to the small perturbations caused by the brown dwarf 
Bb which has the mass of 18.64 Mj (Tab. 4) and consid- 
ering that the Hill stability function is not general enough 
to describe this particular system. Thus it is important to 
note that the establishment of stability using this method 
will be restricted to three body systems. Since significant 
results from this calculation were not produced, we stopped 
the simulation at 8 million years. Thus we became moti- 
vated to modify the HD 41004 system into an ERTBP (Sect 
3.2) and measured the HS values for the planet Ab, shown 
in Fig. 5c. The HS value is constant, positive and demon- 
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Figure 3. Variation of Lyapunov time series for the giant planet as a 4 body problem in (a) HD 41004 Ab and (b) HD 41004 Bb. 
Variation of Lyapunov time series for the giant planet Ab as a 3 body problem (c). Each figure was simulated for 1 x 10^ years. 



strates quasi-periodic oscillations every three million years. 
The amplitude of oscillations reaches a maximum every six 
million years. 

The planet in the HD 41004 system displayed a positive 
stability for the coplanar case irrespective of our choice of 
model. 



3.4 Analysis of MEGNO maps 

The MEGNO indicators are generated over a 150,000 initial 
conditions in eccentricities and semi-major axes for the re- 
spective planets within the selected binaries. In Figs. 6, 7, 8 
and 9 MEGNO maps for different binary eccentricity were 
simulated for 100,000 years. The cross hair in each subplot 
represents the osculating orbit of planet Ab for the respec- 
tive binary (see tables 1 to 4) . The colour bar on top of each 
map indicates the strength in the value of MEGNO (< Y >). 
The blue colour denotes regions of quasi-periodicity and the 
yellow indicates regions of chaos. 

For different eccentricity values of the binary in 7 
Cephei, Fig. 6, the MEGNO indicator shows a clear dis- 
tinction between quasi-periodic and chaotic regions. Within 
the observational value, etin = 0.4085 (Fig. 6a), the planet 
unmistakably demonstrates stable orbits. When the binary 
eccentricity, ebin=0.2 (Fig. 6b), is decreased the cross hair 
is completely inside the quasi-periodic region, hence increas- 
ing the orbital stability. Conversely, as the eccentricity of the 



binary orbit is increased the location of the chaotic mean- 
motion resonances (yellow spikes at constant semi-major 
axis) are shiftet to lower semi-major axis of the planet (Fig. 
6c). 

Figure 7a shows a stable region for the planet in Gliese 
86 system. When the etin is reduced to 0.2 the quasi peri- 
odic region increases (Fig. 7b), indicating that the system 
gets more stable for less eccentric orbits. In Figs. 7b and 
c, we have only considered the semi-major axis interval to 
be [1,4] AU because the observed region [0,1] AU remains 
unchanged with different choices of eccentricity. As a result 
this causes the cross hair to disappear and displays a more 
chaotic region to the right. For a circular binary orbit almost 
all the (a,e) space is quasi-periodic/regular. This is also true 
for the inner region where we have the planet. When the etin 
is increased to 0.6, the chaotic yellow region grows but with- 
out endangering the orbital stability of the observed planet. 

Figure 8a shows a locally stable region for the planet 
in HD 196885 system. The cross hair in the map is located 
right at the edge of the stable region. A small change in 
semi-major axis of the planet can divert the planet towards 
the chaotic region loosing global stability. Quasi-periodicity 
increases with a decrease in the binary eccentricity. Fig. 8b, 
while the system becomes chaotic with increment in binary 
eccentricity. Fig. 8c. 

The planet (Ab) in HD 41004 binary system. Fig. 9, 
has stable orbits even when the binary eccentricity varies 
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Figure 4. Variation of Hill Stability for the giant planet in (a) 7 Cephei, (b) Gliese 86 and (c) HD 196885 simulated for 1 x 10^ yrs. 
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Figure 5. Variation of Hill Stability for the giant planets in HD 41004 as 4 body system (a) and (b), and modified 3 body system (c) 
simulated for 1 x 10^ yrs. 



from to 0.6. Although the planet has a relatively high 
eccentricity (0.39), it resides within the quasi-periodic region 
and maintains global stability. 

3.5 Comparison of Indicators 

The maximum Lyapunov exponent for the planet in 7 
Cephei demonstrates stable orbits. The positive value and 
non-decreasing global trend of Hill stability time series pro- 
vides the necessary evidence. The HS times series also indi- 
cates that the system is quasi-periodic every 5 million years. 
Since the calculation for MLE is limited to 1 million years, 
we are unable to notice any periodicity in the spectra as seen 
in the HS time series but it does also indicate trends toward 
stability. The cross hair in Fig. 6a lies in the quasi periodic 
region of the MEGNO map. Thus, 7 Cephei binary system 
is stable and supplements two of our earlier results. 

The MLE time series. Fig. 2b, calculated for the planet 
in Gliese 86 system shows an oscillating behaviour every 0.3 
million years. The HS time series. Fig. 4a, demonstrates a 
similar oscillation pattern and confirms the comparability of 
the two indicators. The cross hair seen at the periodic re- 
gion of MEGNO map (Fig. 7a) again confirms the stability; 
however, no direct information could be extracted regarding 
the oscillation. 

The MLE, HS time series, and the MEGNO map for 
the planet in HD 196885 system indicate stable orbits. Figs. 
2c, 4c, 8a respectively. The HS spectra is quasi-periodic and 



oscillates every 2 million years. This cannot be deduced from 
the MLE spectra. The MEGNO map on the other hand 
displays a quasi-periodic nature right at the nominal value 
of semi- major axis 2.6 AU. 

The constantly decreasing MLE time series, the posi- 
tive HS time series, and the MEGNO map cross hair is a 
clear indication of stability of the planet in the modified HD 
41004 system (Figs. 3c, 5c, 9a). The time series obtained 
from the MLE and HS values for 1 million years both dis- 
play small oscillations. The reduction of the system from 4 
body to 3 body problem does not affect our MLE calcula- 
tions but significantly affects the HS function and MEGNO 
indicator. The HS values of the planets in the 4 body are 
noisy and produce unusable results while in the reduced sys- 
tem we obtain the oscillating HS values every 3 million years 
which demonstrates the quasi-periodicity of the planet. The 
MEGNO maps for the reduced system indicate the clear 
stability of the planet. 

3.6 Runtime Comparison 

Runtime comparison is vital while comparing the efficiency 
of various indicators due to the fact that timing is crucial 
part of numerical integration. The CPU time depends on the 
computing architecture and integration scheme. We used a 
Yoshida sixth order symplectic integration scheme to calcu- 
late the MLE time series and a Gragg-Bulirsch-Stoer inte- 
gration scheme from Grazier et al. (1996) and Hairer et al. 
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(a) etin = 0.4085 (b) = 0.2 (c) etin = 0.6 

Figure 6. MEGNO Maps for the planets in 7 Cephei simulated for 100,000 years. 
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Figure 7. MEGNO Maps for the planets in Gliese 86 simulated for 100,000 years. 



Runtime 


MLE 


MEGNO 


Hill Stability 


1000 years 


00:13:27 


00:00:01 


00:03:55 


10,000 years 


02:16:15 


00:00:04 


00:35:30 


1,000,000 years 


~ 22 Days 


00:05:25 


~ 3.5 Days 



Table 5. Runtime comparison for different chaos indicators in 
time format of HH:MM:SS. The MLE, MEGNO and HS are com- 
puted for one initial condition for 7 Cephei (Table 1) over given 
period of time and the net CPU run times are compared. 



(1993) to calculate the HS time series and MEGNO maps 
respectively. We believe that one of the reasons for the run- 
time variation is due to the different implementations of the 
Gragg-Bulirsch-Stoer integration scheme. Our runtime com- 
parison is made for 7 Cephei with the nominal initial con- 
ditions given in Table 1. The CPU time to execute each of 
the indicators (MLE, HS time series and MEGNO maps) 
for three different time periods are displayed in Table 5. An 
extensive comparison of runtime is beyond the scope of this 
paper and subject to future studies. 

The Maximum Lyapunov exponent and the Hill stabil- 
ity time series were calculated using Intel(R) Xenon(R) CPU 
E5345 @ 2.33GHz and the MEGNO maps were generated us- 
ing a similar computing architecture, Intel(R) CPU(X5355) 
@ 2.66GHz. The CPU run times shown in Table 5 indi- 
cates that it requires longest amount of time to generate 
MLE time series. The HS time series takes the modest time. 



MEGNO happens to be the fastest to calculate compared 
to two others for a given initial condition (IC). The results 
discussed in this paper were obtained by using one IC for 
MLE and HS time series and 150,000 ICs for MEGNO maps. 
For what we want to achieve from the HS time series, single 
IC would suffice which could be obtained in relatively short 
amount of time. 



3.7 Implications on Habitability 

Based on our previous discussion we know that each planet 
in the four binary stellar systems, 7 Cephei, Gliese 86, HD 
196885 and HD 41004 exists within a stable orbital config- 
uration. Therefore, the existence of other Earth-like planets 
in these systems is a plausible assumption. Most of these 
binary systems have only one planet detected around the 
primary star. Our stability calculations using the MEGNO 
indicator shows a probable region where other possible plan- 
ets could exist. 

For the planet (Ab) in 7 Cephei, eccentricity-semimajor 
axes cross hair lies at [0.115, 1.94] AU. The habitable zone 
(HZ) of the primary star is approximately between [3.05-3.7] 
AU (Haghighipour 2006). The quasi-periodic regions in the 
MEGNO map. Fig. 6a, clearly shows some stability islands 
between [3.0-3.5] AU. For other binary systems, Gliese 86, 
HD 196885 and HD 41004, the HZ is approximated between 
[0.018-1.214] AU, [1.333-2.660] AU and [0.762-1.497] AU, re- 
spectively (Jones et al. 2006). Our MEGNO maps for Gliese 
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Figure 8. MEGNO Maps for the planets in HD 196885 simulated for 100,000 years. 
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Figure 9. MEGNO Maps for the planets in HD 41004 simulated for 100,000 years. 



86 and HD 41004 (Figs. 7a and 9a) display a continuous HZ 
well inside the limits mentioned above. The MEGNO map 
for HD 196885 (Fig. 8a) also displays stability islands within 
the HZ. An Earth-mass test particle can be placed in this 
region and checked for the stability of the system, and even- 
tually looking for the planets in the continuous habitable 
zone. 



4 CONCLUSIONS 

We have applied various chaos indicator techniques in order 
to study the stability of S-type planets in the binaries that 
are less than 25 AU apart. With the time series obtained 
from the maximum Lyapunov exponent and the Hill stability 
function and the maps from the MEGNO indicator, we have 
shown that all the systems are in a stable configuration. 

Our MLE times series for all four planets converge expo- 
nentially which resembles stable planetary orbits. With the 
MLE oscillating every 0.3/0.2 million years in the case of 
Gliese 86 Ab/HD 196885 Ab could be indicative of the near 
resonance behaviour. Similar behaviour is observed in the 
MLE time series of HD 41004 Ab when reduced to 3 body 
problem. Our quest of stability analysis using the MLE time 
series is unaffected by choice of different models of the bi- 
nary system, HD 41004. 

The MEGNO chaos indicator has been effective in 
determining the quasi-periodic regions. The location of 
eccentricity-semimajor axis cross hairs in the MEGNO maps 



for the planets in Gliese 86 and HD 41004 (3 body) systems 
(Figs. 7a and 9a) are well inside the blue region which resem- 
bles the global stability of the planets. The orbital stability 
of these planets are also unaffected by the highly eccentric 
(ebin = 0.6) binary orbits. For the 7 Cephei and HD 196885 
systems (Figs. 6a and 6b), the cross hairs are located in the 
teeth between the chaotic and quasi periodic regions. This 
resembles the local stability of the planets. These planets do 
not survive if the binary orbits are highly eccentric. 

The chaotic and quasi periodic regions observed in the 
eccentricity vs. semi-major axis grid of MEGNO maps were 
the advantage over our two indicators, MLE and HS time 
series, and we used it to explain the possible existence of 
planets in the habitable zone. Each of our chosen binary 
systems shows a possibility of a planet residing within the 
respective habitable zone. 

The Hill stability time series for a planet is successfully 
measured using the potential obtained from numerical inte- 
gration of orbital parameters in elliptic restricted three body 
problem. The measure of the Hill stability is highest at 2.5 
million years and 3 million years for the case of 7 Cephei, 
and this feature is observed to repeat every 5 million years. 
The HS values, nevertheless, stay positive. Similar is the 
case for HD 196885 and HD 41004 (3 body) systems where 
the maximum HS is measured every 2 million years and 3 
million years respectively. The HS values for the planet in 
Gliese 86 system first rises up to 4 and changes very little 
in time. This planet is less perturbed by the secondary star. 
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The analysis of MEGNO indicator for Gliese 86 Ab shows 
a similar results where the planet demonstrates global sta- 
bility. These stability values, for any given system, can be 
extrapolated beyond our 10 million years integration time 
and the stability criteria would still hold true. The HS val- 
ues of 7 Cephei Ab, Ghese 86 Ab, and HD 41004 Ab agree 
with the calculations made by Szenkovits & Mako (2008) 
when the orbital inclination is zero for each system. 

We are able to successfully test the efficacy of Hill sta- 
bility against LEs and MEGNO chaos indicators. Direct 
comparison of stability shows that the Hill stability can be 
set as one of three stringent criterion in the study of sta- 
ble/unstable nature of an planetary orbit in the ERTBP. 
Numerical simulations for the HS function is faster than the 
MLE but slower than the MEGNO indicator. Our results 
show that the HS indicator is comparable to the other in- 
dicators and through our qualitative efficiency analysis we 
have demonstrated that HS can be used as an alternative to 
MLE. 
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